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Three improvements in reduction and com- 
putation of elliptic integrals are made. 1 . 
Reduction formulas, used to express many 
elliptic integrals in terms of a few stan- 
dard integrals, are simplified by modifying 
the definition of intermediate "basic inte- 
grals." 2. A faster than quadratically con- 
vergent series is given for numerical 
computation of the complete symmetric el- 
liptic integral of the third kind. 3. A se- 
ries expansion of an elliptic or hyperelliptic 
integral in elementary symmetric func- 
tions is given, illustrated with numerical co- 
efficients for terms through degree seven 



for the symmetric elliptic integral of the 
first kind. Its usefulness for elliptic inte- 
grals, in particular, is important. 
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Foreword 



Elliptic integrals have many applications, for example 
in mathematics and physics: 

• arclengths of plane curves (ellipse, hyperbola, 
Bernoulli's lemniscate) 

• surface area of an ellipsoid in 3-dimensional space 

• electric and magnetic fields associated with ellipsoids 

• periodicity of anharmonic oscillators 

• mutual inductance of coaxial circles 

• age of the universe in a Friedmann model 

These applications are mentioned in the chapter on el- 
liptic integrals, written by B. C. Carlson, that will ap- 
pear in the NIST Digital Library of Mathematical Func- 
tions . 

The DLMF is scheduled to begin service in 2003 
from a NIST Web site. A hardcover book will be pub- 
lished also. These resources will provide a complete 
guide to the higher mathematical functions for use by 
experienced scientific professionals. The book will 
provide mathematical formulas, references to proofs. 



references to extensions and generalizations, graphs, 
brief descriptions of computational methods, a survey of 
useful published tables, and sample applications. The 
Web site will include, in addition, interactive visualiza- 
tions of 3-dimensional surfaces, a mathematics-aware 
search engine, a downloading capability for equations, 
live links to Web sites that provide mathematical soft- 
ware, and a limited facility for generating tables on 
demand. 

The DLMF is modeled after the Handbook of Mathe- 
matical Functions, published in 1964 by the National 
Bureau of Standards with M. Abramowitz and I. A. 
Stegun as editors. This handbook has been enormously 
successful: it has sold more than 500,000 copies, its 
sales remain high, and it is very frequently cited in 
journal articles in physics and many other fields. But 
new properties of the higher functions have been devel- 
oped, and new functions have risen in importance in 
applications, since the publication of the Abramowitz 
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and Stegun handbook. More than half of the old hand- 
book was devoted to tables, now made obsolete by the 
revolutionary improvements since 1964 in computers 
and software. The need for a modern reference is being 
filled by more than 30 expert authors who are working 
under contract to NIST and supervised by four NIST 
editors. The writing is being edited carefully to assure 
consistent style and level of content. 

Elliptic integrals have long been associated with the 
name of Legendre. Legendre's incomplete elliptic inte- 
grals are 



which follows is an example. It is a further development 
of material that appears in Carlson's DLMF chapter on 
elliptic integrals. 

Daniel W. Lozier 

NIST Mathematical and Computational Sciences 

Division 



1. Simplified Formulas for Reducing 
Elliptic Integrals 



F{cl>,k) = 



d0 



vr 



k^ sin" e 
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E{(l),k)=\ Vl - fe'sin'ede, 
Jo 



and 

n{cf),a\k): 



dO 



Vl - fe'sin'0 (1 - a'- sin' 6)' 



The complete forms of these integrals are obtained by 
setting (j> = 17/2. 

Over a period of more than 35 years, Carlson has 
published a series of papers that provide valuable new 
mathematical and computational foundations for the 
subject in terms of the symmetric elliptic integrals 



RF{x,y,z)- 



f 



df 



A large class of elliptic integrals can be written in the 
form 



J V i^ I i^ 1 



1) 



where m = (mi, . . ., m„) is an « -tuple of integers (posi- 
tive, negative or zero), x>y, h = 3 or 4, n^h, and the 
different linear factors are not proportional. The a 's and 
b 's may be complex (with the b 's not equal to zero), but 
the integral is assumed to be well defined, possibly as a 
Cauchy principal value. In particular the line segment 
with endpoints a, + bix and a, + biy is assumed to lie in 
the cut plane C\(-=c, 0) for 1 < / < /;. 

We write m = S^Li mjCj, where e^ is an « -tuple with 1 
in the jth position and O's elsewhere, and we define 
= (0, . . ., 0). Reference [1] gives a general method of 
reducing /(m) to a linear combination of "basic inte- 
grals," defined as 7(0), /(— Cy), 1 <7 < n, and (if h = 4) 
I(ek), 1 < /: < 4. A simple example of such a reduction 



RD(x,y,z) = - 



dt 



2j,, s{t)it + zy 



bi /(e, - e,) = dj, /(-e,) + bj 1(0), 



i,J- 



[1, 



(1.2) 



where 



Rj(x,y,z,p)-- 



dt 



s(t){t + p) 



s(t) = \/t + x \/t + y 



The complete forms are obtained by setting x = 0. In 
comparison with Legendre's integrals, Carlson's inte- 
grals simplify the reduction of general elliptic integrals 
to standard forms and open the way to efficient compu- 
tations by application of a duplication theorem. 

One of the purposes of the DLMF project is to stim- 
ulate research into the theory, computation and applica- 
tion of the higher mathematical functions. The paper 



where 4, = (ijbi — aibj . This equation reduces all 36-1-72 
integrals in Ref. [2], Eqs. (3.142) and (3.168) and also 
the 18 integrals in Ref. [2], Eq. (3.159) after taking x^ 
as a new variable of integration. The basic integrals are 
expressed in terms of symmetric standard integrals Rf, 
Rd, and 7?, in Ref. [1], Sec. 4. 

The general method first reduces /(m) by Ref. [1], 
Eq. (2.19) to integrals in which m has at most one 
nonzero component and then uses two recurrence rela- 
tions Ref. [1] Eqs. (3.5) and (3.11) for further reduction 
to basic integrals. For example, the simplest recurrence 
relation is Ref. [1], Eq. (3.11): 



fo,?/(^e,) = 2rK^r/(re,), q&N. (1.3) 
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The b 's appear also in the other two formulas and there- 
fore in all reduction formulas, sometimes in great profu- 
sion if m is considerably more complicated than in Eq. 
(1.2). 

It has been found that the b 's disappear from all three 
formulas, and therefore from all reduction formulas, if 
we define 



/(m) = I{m)IB , A (m) = A (m)/5 , 

Here A (m) is the algebraic function 
A(m)=/„(x) -fAy), 



(1.4) 



(1.5) 



where /m(f) is the integrand of Eq. (1.1). Note that 
/(0)=/(0). 

For example, Eq. (1.2) becomes 



/(e,-e,) = r,,/(-e,)+/(0), 
and Eq. (1.3) becomes 



(1.6) 

hq^i) = l.[\]rri{r^i\ ? G N. (1.7) 

The other recurrence relation, Ref. [1], Eq. (3.5), be- 
comes 

h 

2(2^ + r)Ei,-,{rij, . . ., nj)Iiiq - 1 + r)ej) 
= 2A(qej + ^t), q&Z, l<;<n (1.8) 



where Ei,-r is the elementary symmetric function de- 
fined by 



[1(1 + tr.j) = ^t%iry, . . ., r.j), 



(1.9) 



whence Eh = Oifl^j^h because r,j = 0. 

The remaining formula, Ref. [1], Eq. (2.19), becomes 

M n -rrtj 

/(m) = J,CM-,ik)I{q^,) + J,A^C„„,,(i)I{-qtd, 

(1.10) 

where M = 2]Li rtij and each sum is empty if its upper 
limit is less than its lower limit. The first term on the 
right is independent of k, which is usually best chosen 
so that 1 < k ^ h. The coefficients are defined by 



A.=n,-, co(o=i, c..(/)=E "^^'i::::f'\ 






tti ! • • • a, ! 



V±p(i) = ^im,rr, p^l, (1.11) 



where upper (lower) signs go together and the first sum 
extends over all nonnegative integers ai, . . ., a, such that 
ai + 2a2 + ... + sus = s. A recurrence relation for the 
coefficients is 



^C±,(/) = 2/'^±,(/)C±(,_„)(/), .v>l. (1.12) 

p=i 



2. Algorithms for Complete Elliptic 
Integrals of the Third Kind 

Complications formerly encountered in numerical 
computation of Legendre's complete elliptic integral of 
the third kind were avoided by defining and tabulating 
Heuman's Lambda function (for circular cases) and a 
modification of Jacobi's Zeta function (for hyperbolic 
cases). For example, the method of Ref. [3] was later 
superseded by Bartky's transformation and its applica- 
tion by Bulirsch [4] to his complete integral 
cel{kc, p, a, b). Bartky's transformation for the com- 
plete case of the symmetric integral of the third kind, 
obtained from Ref. [5], Eq. (4.14) with the help of 

(37i/4)Rdy, z,p) = 3Rf{0,y, z) - pRAO,y, z,p). 



(Ti/2)RK(y,z) = RF(0,y,z), 
can be written as 



(2.1) 



Rj(0, gl, al p^) = s„Rj(0, glu alu pL) 

+ (3/2p^)Rf(0, gLu al,), (2.2) 

where a„ , g„, Pn are positive, s„ = (pH — af,gl )l%Pn , and 



(1,1+1 — ' 



an + g,, 

2 ' 



4-1 = va„g„, p„+i = 



N. 



Pn + Ongn 
'2.Pn 



(2.3) 



As n ^ 00, a„ and g„ converge quadratically to 
Gauss's arithmetic-geometric mean, M(flo, go), and 

7?f (0, gf,, al) = Ti/lMiao, go), nGN, (2.4) 

by Ref. [6], Eqs. (6.10-8) and Eq. (2.1). It follows from 



415 



Volume 107, Number 5, September-October 2002 

Journal of Research of the National Institute of Standards and Technology 



(2.3) that 



P„+l - g„+l = (Pn- g„+lf/2p„. 



(2.5) 



Since g„+i converges quadratically to M(flo, go), so does 
p„ , and s„ converges quadratically to 0. Iteration of Eq. 
(2.2) gives 

RAO, gl al pi) = ^ RjiO, gf, , ai , p^ ) 



Po 
^ 4piMiao, go) io 



Vk'here 



(2.6) 



where we have used Eq. (2.4) and chosen x, = ai in Ref. 
[7], Eq. (4.6). Substitution of Eq. (2.9) gives 



Rj(0, gl al -qi) = 



-3tt 



4M(flo, go){qS + flo) 



2 + 



flp - go 

?o + go to 



2g»- 



(2.12) 



Equation (2.3) and the second line of Eq. (2.9) still 
apply, with po given by Eq. (2.11). 

For the complete case of Legendre's third integral, 

U(a\ k) = (ay3)Rj(0, k'\ 1, 1 - a') + Rf{Q, k'\ 1), 

(2.13) 



20=1, -y = -2 , m>l, (2.7) 

Po Pm 



and therefore 

Q„+i _ .?„/?,? _ 1 p„^ - fl„g„ 
Qn pL\ 2 /7„^ + a„g„ ' 



(2.8) 



Letting n ^ co in Eq. (2.6), we find, for positive 

flo, go, Po, 



where k'^ = I — k^, Eq. (2.9) implies 



nia\k) 



2.T^2C. 



4M(l,k') 
CO <^^ < 1, -00 < a^ < 1 



(2.14) 



where Eq. (2.3) applies with ao=l, go = k', and 
Po' = 1 - a\ 

If a' > 1, Eq. (2.12) gives the Cauchy principal value, 



RAO, gl ai, pi) ■■ 



3tt 



4piM(ao, go) t:!) 



2 a. 



^(a^yfe): 



-TTr 



4M(l,fe')(a - 'tot; 



2 a, 



Go=l, G„.4g„.„, ^. = fj^, (2.9) 



where e„ converges to quadratically and Q„ converges 
to faster than quadratically. 

If Pa = flo, Eq. (2.9) becomes 



-00 <fc^ < 1, 1 < a' < oo^ 



(2.15) 



RD(0,go,ao) = -7-2 



3-77 



2 a, 



4floM(flo, go) ^ 



Go=l, G«..=^ae„, e„ = ^--f=, (2.10) 

^ "« + gn 



where /?d is a complete integral of the second kind, 
symmetric in only its first two arguments. If < flo — go 
then - 1 < Co < 0, but e„ > and g„ < for n > 1. 

If the last variable of R, is negative, the Cauchy prin- 
cipal value is given by 

(qi + ai)Rj{0, gi, ai, -qi) = (pi - ai)RAO, gi, ai, pi) 

3tt 



2M{ao, go) ' 



pi = aiiqi + gi)l(qi + ai), (2.11) 



where Eq. (2.3) and the second line of Eq. (2.9) apply 
with ao = 1, go = k', and po = 1 — k^/a^. 



3. Expansion in Elementary Symmetric 
Functions 

The duplication method of computing the symmetric 
elliptic integrals Rp and R, (including their degenerate 
cases Re and Rd) consists in iterating their duplication 
theorems until their variables are nearly equal and then 
expanding in a series of elementary symmetric functions 
of the small differences between the variables. In the 
absence of a duplication theorem the method is useful 
for a hyperelliptic integral only if the variables are 
nearly equal. The series is truncated to a polynomial of 
fixed degree; the higher the degree, the fewer duplica- 
tions are needed for a desired accuracy of the result but 
the larger the number of terms to be calculated. No tests 
have been made to determine an optimal compromise, 
which would depend in large part on the speed of 
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extracting square roots in the duplication theorem and 
therefore on the equipment used. In Ref. [8] polynomi- 
als of degree five were chosen for simplicity, but later it 
seemed worthwhile to increase the degree for the com- 
paratively slow computation of R,, and Ref. [9], Eq. 
(A. 11) gives the terms of degree six and seven for R,. 
The change in speed would be significant only if a very 
large number of computations were performed, since 
the result of a single computation is returned with no 
delay apparent to the eye. We shall give here the corre- 
sponding terms for Rf and the general form of the in- 
finite series of which the polynomials are truncations. 
Any 7? -function R-a{bi, ..., bk', Zi,.-.,Zk) (see Ref. 
[6]) in which all the Z? -parameters are positive integral 
multiples of a single number j3 can be rewritten with 
repeated variables and all b 's equal to j3 . An example is 

Rj(x , y , z, p) = R-ij2(h h h 1; x,y,z,p) 

= R-y2(iiii-2;x,y,z,p,p), (3.1) 

and RD{x,y,z) is the case with p = z. Therefore we 
consider 

7?-„(;8,..., )S; zi, ...,z„) 
1 



;fr--'fl 

•'0 1=1 



(t + z.r^dt 



B(a, nf3 — a ) . 

= A-''R-,{I3,..., I3; zi/A,..;Z„/A), (3.2) 

where B is the beta function and we have used the 
homogeneity of R to divide the variables by their arith- 
metic average, 

A=-Sz,.. (3.3) 

" 1=1 

The relative difference between A and z, is 



Zi = - 



= 1 



Zj 

A' 



(3.4) 



whence 



R-,.{l3,...,fi;zu-.;Zn) 
■■A-''R^,(I3,...,I3; 1 -Zi, ..., 1 - Z„). 



(3.5) 



Because the function is symmetric in the Z's, it can 
be expanded in elementary symmetric functions 
E,„ = E„iZx, . . ., Z„ ) defined by 



Y\i\+tZd = J,fE, 



(3.6) 



Applying Ref. [10], Eqs. (A.5) and (A.12) to the right 
side of Eq. (3.5), we find 



R-Aii,...,ii\zx 



.Zn) 



=A-E 



{a)N 
{nji), 



T^il3,..., 13 ;!,,... ,Z„), 



(a). 



A-"27Sf2(-lfn^). 



(3.7) 



(3.8) 



where (fl)^ is Pochhammer's shifted factorial, M = SJLi 
OT,, and the inner sum (representing Tn) extends over all 
nonnegative integers OTi, ...,ot„ such that mi + 2m2 
+ ■■• + nm„ = N. Reference [10], Eq. (A. 6) provides the 
recurrence relation 



NT^ + 2(- 1)X'-|S + N- r)EJ^^, = 0, 



(3.9) 



where To = 1 and T^-r = if r> N, whence Ti = 0. The 
term with r= 1 is missing because Eqs. (3.3) and (3.4) 
imply El = 0, which greatly simplifies T^. 

Let IZI = max, IZ,I and A = max{lal, 1}. The series 
[Eq. (3.7)] converges absolutely if /3 > and IZI < 1, and 
the truncation error rK resulting from neglect of terms of 
degree N ^ K can be shown to satisfy 



Ir^l 



(lal)glZI*^ 
KHl - IZlf 



(3.10) 



Each application of a duplication theorem decreases by 
a factor of four the differences between the variables, 
therefore the difference A — zi, and ultimately IZI as A 
approaches a limit. It is easy to determine whether an- 
other duplication is needed before using the truncated 
series to achieve the desired accuracy. 

If a = jS = 1/2 and n = 3, the series [Eq. (3.8)] with 
£1=0 can be rearranged as 



Rf(zi,Z2,z,)=A-"'J,J, 



i-mkU EiEl 



■^^ Ar+6s+\r\s\ 



(3.11) 



This double series is convenient for obtaining numerical 
coefficients but requires lEal + l^sl < 1 for absolute con- 
vergence. Keeping together all terms of the same degree 
in the Z's, as in Eq. (3.8), we find 



Rf{zuZ2,z,)=A-"\\ 



10 



14' 



£2 + -r7 £3 + :r7 El 



24' 



44 



E2E, -^El + ^EJ + ^ EiE, + r,), (3.12) 
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where 



4. References 



\n\< 



0.2\Zf 

1 - izr 



IZI = max 



< 1, 



A = (zi+Z2 + Zi)l3. 



(3.13) 



One more duplication would (ultimately) have de- 
creased I rsl by a factor of 4^ = 65 536. 

For convenient reference we restate the corresponding 
result for Rj and Ro [see Eq. (3.1) and the discussions 
preceding it]: 



R-viii 



;z., 



.Zs) 



■A (1 — [4 £2 + 6-^3 + 88 -£'2 ~22 -£'4 



9_ 

52 



E,E, + — E, 



j^Ei + -Ei + -E,E, 



+ ^EiE, 
272 
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E3E4 



68 



E2ES + rs). 
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where 



lrsl< 



3.4izr 

(1 - \Z\f 



IZI = max 



< 1, 



A = (ZiH-...-HZ5)/5. 



(3.15) 



About the author: B. C. Carlson was formerly a pro- 
fessor of physics and is currently a professor emeritus 
of mathematics at Iowa State University and an associ- 
ate at the Ames Laboratory of the U.S. Department of 
Energy. 



The duplication theorems for these two functions are 
more complicated than that for Rp; see Ref. [10] and 
Ref. [9], Appendix. 
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